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ABSTRACT 

To  test  the  hypothesis  that  high  quality  3D  Earth  models  will  produce  seismic  event  locations  which  are  more 
accurate  and  more  precise,  we  are  developing  a  global  3D  P  wave  velocity  model  of  the  Earth’s  crust  and  mantle 
using  seismic  tomography.  In  this  paper,  we  present  the  most  recent  version  of  our  model,  SALSA3D  (SAndia  LoS 
Alamos)  version  1 .4,  and  demonstrate  its  ability  to  reduce  mislocations  for  a  large  set  of  realizations  derived  from  a 
carefully  chosen  set  of  globally-distributed  ground  truth  events. 

Our  model  is  derived  from  the  latest  version  of  the  Ground  Truth  (GT)  catalog  of  P  and  Pn  travel  time  picks 
assembled  by  Los  Alamos  National  Laboratory.  To  prevent  over- weighting  due  to  ray  path  redundancy  and  to 
reduce  the  computational  burden,  we  cluster  rays  to  produce  representative  rays.  Reduction  in  the  total  number  of 
ray  paths  is  >  55%.  The  model  is  represented  using  the  triangular  tessellation  system  described  by  Ballard  et  al. 
(2009),  which  incorporates  variable  resolution  in  both  the  geographic  and  radial  dimensions.  For  our  starting  model, 
we  use  a  simplified  two  layer  crustal  model  derived  from  the  Crust  2.0  model  over  a  uniform  AK135  mantle. 
Sufficient  damping  is  used  to  reduce  velocity  adjustments  so  that  ray  path  changes  between  iterations  are  small.  We 
obtain  proper  model  smoothness  by  using  progressive  grid  refinement,  refining  the  grid  only  around  areas  with 
significant  velocity  changes  from  the  starting  model.  At  each  grid  refinement  level  except  the  last  one  we  limit  the 
number  of  iterations  to  prevent  convergence  thereby  preserving  aspects  of  broad  features  resolved  at  coarser 
resolutions.  Our  approach  produces  a  smooth,  multi-resolution  model  with  node  density  appropriate  to  both  ray 
coverage  and  the  velocity  gradients  required  by  the  data.  This  scheme  is  computationally  expensive,  so  we  use  a 
distributed  computing  framework  based  on  the  Java  Parallel  Processing  Framework,  providing  us  with  ~400 
processors.  Resolution  of  our  model  is  assessed  using  a  variation  of  the  standard  checkerboard  method,  as  well  as  by 
directly  estimating  the  diagonal  of  the  model  resolution  matrix  based  on  the  technique  developed  by  Bekas,  et  al. 

We  compare  the  travel-time  prediction  and  location  capabilities  of  this  model  over  standard  ID  models.  We  perform 
location  tests  on  a  global,  geographically-distributed  event  set  with  ground  truth  levels  of  5  km  or  better.  These 
events  generally  possess  hundreds  of  Pn  and  P  phases  from  which  we  can  generate  different  realizations  of  station 
distributions,  yielding  a  range  of  azimuthal  coverage  and  proportions  of  teleseismic  to  regional  arrivals,  with  which 
we  test  the  robustness  and  quality  of  relocation.  The  SALSA3D  model  reduces  mislocation  over  standard  ID  akl35, 
especially  with  increasing  azimuthal  gap.  The  3D  model  appears  to  perform  better  for  locations  based  solely  or 
dominantly  on  regional  arrivals,  which  is  not  unexpected  given  that  akl35  represents  a  global  average  and  cannot 
therefore  capture  local  and  regional  variations. 
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INTRODUCTION 

The  evolving  U.S.  monitoring  needs  require  accurate  location  of  smaller  and  smaller  events.  Given  the  limited 
station  coverage  of  a  typical  global  monitoring  network,  such  as  the  International  Monitoring  System  (IMS),  it  is 
likely  that  such  events  will  be  detected  by  a  small  number  of  stations,  quite  possibly  with  a  poor  network  geometry. 
For  locating  such  events,  it  is  essential  to  have  extremely  high  fidelity  travel  time  predictions,  particularly  at 
regional  distances  where  lateral  heterogeneity  can  be  significant.  Accurately  accounting  for  lateral  heterogeneity 
implies  using  3D  models  of  the  Earth  to  calculate  travel  times,  but  relatively  few  of  the  available  3D  Earth  models 
are  appropriate  for  high  fidelity  travel  time  prediction,  and  it  is  unclear  whether  any  of  those  that  are  appropriate 
actually  improve  event  locations.  Thus,  we  are  developing  our  own  global  3D  P  wave  velocity  model  of  the  Earth’s 
crust  and  mantle  —  SALSA3D  (SAndia  LoS  Alamos  3D)  —  using  seismic  tomography  based  on  a  carefully 
assembled  data  set  of  P  phase  travel  times  collected  by  LANL  over  the  past  decade.  An  important  difference 
between  our  effort  and  previous  efforts  is  that  our  model  was  produced  specifically  for  improving  event  location;  the 
model  provides  potentially  valuable  information  about  the  structure  of  the  Earth,  but  this  is  not  our  focus.  Hence,  all 
decisions  about  data  processing  and  tomography  were  made  with  this  goal  in  mind. 

In  this  paper,  we  present  the  most  recent  version  of  our  model,  SALSA3D  version  1 .4,  and  demonstrate  its  ability  to 
reduce  mislocations  for  a  large  set  of  realizations  derived  from  a  carefully  chosen  set  of  globally-distributed  ground 
truth  events. 

DATA  SET 

The  data  used  for  the  tomographic  inversion  was  collected  by  LANL  over  the  past  decade  and  includes  events 
categorized  with  a  ground  truth  (GT)  level  of  25  km  or  better  (Bondar  et  al.,  2004).  Arrivals  for  these  GT  events 
were  merged  from  various  sources  in  order  to  produce  a  single  set  of  arrivals  for  each  event  with  redundancies 
removed,  maximizing  the  available  arrivals  for  each  event  (Begnaud,  2005).  Regional  Pn  and  teleseismic  P  arrivals 
were  extracted  for  use  in  the  tomography  and  resulted  in  ~1 19K  events,  ~12K  stations,  and  almost  12  million 
individual  ray  paths  (Figure  1).  Geographically  focused  component  data  sets  include  USarray  and  the  Deep  Seismic 
Soundings  in  the  former  Soviet  Union  (e.g.,  Li  and  Mooney,  1998).  To  prevent  over- weighting  due  to  ray  path 
redundancy  and  to  reduce  the  computational  burden,  we  cluster  rays  to  produce  representative  rays,  thereby  reducing 
the  total  number  of  ray  paths  by  55%. 


Figure  1.  Stations  and  events  used  for  the  global  tomography  inversion.  These  data  were  collected  by  LANL 
over  the  past  decade  and  are  all  categorized  with  a  ground  truth  level  of  25  km  or  better. 
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TOMOGRAPHIC  INVERSION 

We  use  the  LSQR  algorithm  of  Paige  and  Saunders  (1982)  to  invert  the  data  to  find  the  P  wave  velocity  distribution 
in  the  mantle  that  optimally  reduces  the  misfit  between  observed  and  predicted  travel  times  for  P  and  Pn  phase 
arrivals.  Our  starting  model  for  the  inversion  consists  of  a  modified  version  of  the  Crust  2.0  model  (Laske  and 
Masters,  1997,  Laske  et  al.,  http://igppweb.ucsd.edu/~gabi/crust2.html)  overlaying  an  AK135  mantle.  We  modify 
Crust  2.0  by  combining  their  six  layers  into  just  two.  The  crust  is  damped  heavily  in  the  inversion  so  it  changes  very 
little.  We  retrace  the  rays  at  each  of  16  iterations  using  the  pseudo-bending  algorithm  (Um  and  Thurber,  1987,  Zhao 
and  Lei,  2004,  Ballard  et  al.,  2009).  During  early  iterations,  the  mantle  is  moderately  damped  to  discourage  the 
inversion  from  progressing  too  rapidly  but  then  damping  is  relaxed  in  later  iterations.  Adding  regularization 
constraints  to  the  data  kemal  matrix  is  not  required  during  the  inversion  since  it  is  accomplished  by  beginning  with  a 
coarse  grid  that  is  progressively  refined  several  times  during  the  inversion  (Simmons  et  al.,  2009). 

ADAPTIVE  GRID  REFINEMENT 

P  wave  velocity  in  our  model  is  stored  on  a  3D  grid  of  nodes  that  has  variable  resolution  in  both  geographic  and 
radial  dimensions.  Our  basic  starting  point  is  a  2  dimensional,  multi-level  triangular  tessellation  of  a  unit  sphere 
(Wang  and  Dahlen,  1995,  Ballard  et  al.,  2009).  We  assign  a  tessellation  with  8°  triangles  to  the  lower  mantle,  a 
tessellation  with  4°  triangles  to  the  transition  zone  and  upper  mantle,  and  a  third  tessellation  with  variable  resolution 
to  all  crustal  layers.  The  crustal  tessellation  (not  shown)  has  2°  triangles  in  oceanic  regions,  1°  triangles  in  most 
continental  regions,  and  0.5°  triangles  in  the  United  States  where  the  additional  resolution  is  warranted  due  to  the 
dense  station  coverage  afforded  by  the  US  Array.  The  two  tessellations  in  the  mantle  are  refined,  in  both  the 
geographical  and  radial  dimensions,  during  the  progress  of  the  tomographic  inversion.  Figure  2  illustrates  the  final 
geographic  node  configuration  of  the  tessellation  assigned  to  the  transition  zone  and  upper  mantle.  4°  triangles  in 
relatively  aseismic  parts  of  the  world  remain  unrefined  in  the  final,  model,  but  the  grid  is  refined  down  to  triangles 
as  small  as  0.5°  in  areas  with  dense  seismicity  or  stations,  or  both. 


Figure  2.  The  final  grid  for  the  transition  zone  and  upper  mantle  with  triangle  sizes  ranging  from  0.5°  to  4°. 
The  grid  for  these  layers  in  the  starting  model  consisted  entirely  of  4°  triangles. 

VELOCITY  MODEL 

A  small  portion  of  the  global  range  of  our  model  is  illustrated  in  Figure  3  and  Figure  4.  The  images  reveal  many 
tectonic  features  of  the  Earth.  Major  cratons  in  Siberia,  North  America,  Southern  Africa  and  Australia  are  all 
characterized  by  P  wave  velocities  that  are  relatively  fast  compared  to  AK135.  Tectonically  active  areas  such  as 
western  North  America  and  the  Mediterranean  region  are  relatively  slow,  as  are  the  mid-ocean  ridges  in  the  Pacific, 
Atlantic  and  Indian  Oceans.  In  subduction  zones  beneath  the  Andes,  Sumatra  and  New  Zealand  we  see  fast  oceanic 
slabs  being  subducted  beneath  the  adjacent  continents.  The  stagnant  subducted  slab  beneath  Japan  is  particularly 
evident  in  Figure  4.  Very  fast  material  is  evident  in  the  upper  mantle  beneath  Tibet  and  slow  velocity  anomalies 
appear  in  the  Red  Sea  region. 
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Depth  =  222  km 


Vp  %  Change  from  AK135 
-2-1012 


Figure  3.  Map  of  the  %  difference  in  P  wave  velocity  between  SALSA3D  and  AK135  at  a  depth  of  225  km  in 
the  upper  mantle.  The  black  line  shows  the  position  of  the  approximate  great  circle  path  of  the 
cross  section  shown  in  Figure  4. 
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Figure  4.  Cross  section  through  SALSA3D  at  the  location  shown  by  the  black  line  in  Figure  3. 


RESOLUTION  TESTS 


Images  of  tomographic  models  can  provide  important  insight  into  fundamental  questions  about  the  structure  of  the 
Earth,  such  as  the  fate  of  subducting  slabs.  However,  when  interpreting  these  images,  one  must  consider  the 
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resolution  of  the  tomographic  inversion.  For  large  models  like  ours,  this  can  be  a  computationally  challenging 
problem.  A  standard  way  to  estimate  model  uncertainty  is  to  create  known  detailed  velocity  structures  (often 
’’checkerboard"  patterns),  use  these  to  create  synthetic  observations  for  the  actual  data  paths,  and  then  to  invert  these 
synthetic  observations  to  see  how  much  of  the  known  structure  can  be  recovered.  This  technique  measures  the 
theoretical  resolution  possible  given  the  ray  coverage  (i.e.  it  is  based  entirely  on  ray  geometries).  The  actual 
resolution  is  something  less  than  this  because  the  observations  themselves  will  have  errors. 

We  adopt  the  checkerboard  technique  here,  though  our  pattern  is  slightly  different,  consisting  of  equal  sized 
alternating  velocity  triangles,  so  that  we  can  maintain  the  same  feature  size  regardless  of  latitude  (the  rectangular 
latitude/longitude  checkers  used  by  many  researchers  actually  get  smaller  towards  the  poles).  Our  checkers  represent 
deviations  relative  to  our  final  tomographic  model  so  that  the  ray  paths  will  be  as  realistic  as  possible.  For  the  same 
reason,  our  checkers  have  relatively  small  deviations  (+/-  1%)  relative  to  the  background  model. 

As  expected  given  the  ray  coverage  (Figure  1),  our  model  resolution  varies  strongly  both  laterally  and  with  depth. 
Resolution  is  best  for  Asia,  due  to  the  excellent  coverage  provided  by  the  LANL  GT  data  set.  Resolution  is  best  in 
the  upper  portion  of  the  lower  mantle,  due  to  the  large  number  of  teleseismic  raypaths  in  our  dataset  that  bottom  near 
this  depth. 


Vp  %  Change 

-1  -0.5  0  0.5  1 


Figure  5.  Results  of  the  Chinese  checkerboard  test  at  a  depth  of  1139  km  in  the  shallow  part  of  the  lower 
mantle. 

TRAVEL  TIME  RESIDUAL  REDUCTION 

To  evaluate  the  ability  of  our  model  to  accurately  predict  observed  source-receiver  seismic  travel  times,  we  define  a 
set  of  19  polygons  that  surround  areas  of  contiguous  seismicity  around  the  world  (Figure  6).  For  each  of  the  49 
seismic  stations  in  the  IMS  primary  network  we  compute  the  mean  absolute  residual  (absolute  value  of  observed 
minus  predicted  travel  time)  for  all  of  the  events  in  each  polygon  using  the  AK135  model  (Kennett  et  al.,  1995)  and 
using  SALSA3D.  For  computing  the  AK135  travel  times  we  used  the  Taup  Toolkit  software  package  (Crotwell  et 
ah,  1999).  The  results,  illustrated  in  Figure  7,  indicate  that  SALSA3D  produces  travel  time  residuals  that  are  smaller 
than  those  produced  by  AK135  more  that  92%  of  the  time. 
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Figure  6.  Map  of  the  locations  of  the  49  IMS  primary  network  stations  and  the  19  seismic  regions  selected 
for  travel  time  residual  analysis. 


Figure  7.  Median  absolute  value  of  residual  computed  with  AK135  vs  the  same  value  computed  with 
SALSA3D.  Each  point  represents  a  single  station-region  pair. 

LOCATION  IMPROVEMENT 

To  test  the  location  accuracy  of  the  SALSA3D  model,  we  selected  a  set  of  GT5  or  better  events  from  throughout 
Eurasia  that  contained  enough  Pn  and  P  arrivals  to  be  able  to  create  a  set  of  random  realizations  of  differing  station 
coverages.  42  events  in  Eurasia  were  initially  selected  and  modified  from  the  location  test  set  of  Myers  et  al.  (2010), 
pulling  data  from  the  LANL  GT  database  (Figure  8). 
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For  each  event,  5  random  selections  of  P  and  Pn  arrivals  were  made  for  each  target  number  of  arrivals  from  5  to  50 
(incremented  by  1).  The  total  number  of  realizations  was  9585  (one  event  did  not  have  50  available  phases). 

Our  event  clustering  process  performed  prior  to  the  tomography  significantly  reduces  any  circularity  issues  despite 
the  lack  of  a  formal  event  segregation  procedure. 

We  relocated  each  realization  with  the  standard  ID  AK135  model  (Kennett  et  al,  1995)  and  our  global  SALSA3D 
tomography  model.  For  AK135,  we  used  the  appropriate  ellipticity/station  corrections  and  the  default  model  error. 
For  SALSA3D,  an  estimated  model  error  was  used  based  on  a  summary  standard  deviation  of  residuals  with 
distance.  Results  were  evaluated  based  on  median  mislocations  with  the  total  number  of  phases,  azimuthal  gap,  and 
the  relative  number  of  teleseismic  P  versus  regional  Pn. 

Figure  9  shows  the  median  mislocation  of  the  relocated  random  realizations  with  azimuthal  gap.  At  small  azimuthal 
gaps,  the  median  mislocation  for  the  3D  model  is  slightly  reduced  over  AK135,  but  the  good  station  coverage 
produces  accurate  locations  for  either  model.  As  the  azimuthal  gap  increases,  the  use  of  the  3D  model  results  in  a 
greater  reduction  in  mislocation  over  the  ID  AK135  model  by  as  much  as  a  ~10  km. 

Figure  10  shows  the  median  mislocation  using  different  combinations  of  P  and  Pn  arrivals.  The  maximum  number 
of  any  arrival  set  is  50.  The  SALSA3D  model  has  an  overall  lower  median  mislocation.  The  combination  of  equal 
numbers  of  Pn  and  P  arrivals  results  in  the  highest  median  mislocation  values  for  the  ID  AK135  model.  This 
demonstrates  the  problems  combining  regional  and  teleseismic  arrivals  when  using  a  model  not  designed  for  such  a 
purpose.  The  results  using  the  3D  model  do  not  show  this  higher  median  mislocation  in  the  area  of  equal  P  and  Pn, 
indicating  reduced  problems  with  combining  regional  and  teleseismic  phases.  Notice  that  results  with  fewer  phases 
(<  10)  clearly  demonstrate  improvement  by  the  3D  model,  as  expected. 


Figure  8.  Validation  events  for  location  testing.  42  GT5  or  better  events  were  selected  within  Eurasia  in 
order  to  produce  a  random  data  set  of  varying  station/data  configurations  in  which  to  test  the 
location  accuracy  of  the  3D  model.  The  events  were  selected  from  Myers  et  al.  (2010)  with  a  few 
additional  to  improve  spatial  coverage. 
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Azimuthal  Gap  (degrees) 

Figure  9.  Median  mislocation  versus  azimuthal  gap  for  the  location  improvement  testing  using  the  random 
realizations  of  the  42  events  shown  in  Figure  8. 
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Figure  10.  Median  mislocation  with  varying  number  of  P  and  Pn  arrivals  for  the  random  realization  location 
tests  (bin  size  =2).  The  AK135  model  results  (top  left)  display  higher  mislocation  values  than  the 
SALSA3D  model  (top  right).  The  higher  mislocation  values  at  approximately  equal  numbers  of  P 
and  Pn  arrivals  seen  in  the  AK135  results  are  not  observed  for  the  3D  model.  The  mislocation 
difference  between  the  models  (bottom)  shows  that  the  3D  model  is  performing  significantly  better 
than  AK135  for  the  majority  of  phase  combinations. 
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CONCLUSION 

In  this  paper,  we  have  presented  the  latest  version  of  SALSA3D,  the  Sandia-Los  Alamos  global  tomography  3D  P- 
velocity  model  of  the  crust  and  mantle.  Our  model  is  multi-resolution,  with  the  resolution  determined  dynamically 
as  part  of  the  tomographic  inversion.  Additional  nodes  are  added  where  ray  coverage  is  high  and  velocity  gradients 
are  significant;  other  areas  maintain  the  starting  grid  sampling.  SALSA3D  images  many  well-established  features 
within  the  Earth,  such  as  subducting  slabs  and  continental  roots,  and  is  generally  consistent  with  other  recent  global 
and  regional  models.  Resolution,  as  established  by  checkerboard  tests,  is  highly  variable  both  laterally  and  with 
depth,  with  the  best  resolution  in  Asia  and  Europe  in  the  upper  part  of  the  lower  mantle. 

More  important  than  the  structural  features  in  the  model  is  the  fact  that  SALSA3D  provides  compelling 
improvement  in  travel  time  prediction  and  in  event  location  compared  to  AK135.  Evaluating  performance  for  the 
IMS  primary  network,  travel  time  residuals  grouped  by  station  and  geographic  region  show  clear  improvements  for 
virtually  all  combinations  spanning  a  variety  of  distances  and  tectonic  regions.  Event  location  improvement  is 
evaluated  by  re-locating  multiple  realizations  from  a  set  of  42  well-characterized  events  in  Europe  and  Asia. 
SALSA3D  improves  locations  only  slightly  when  the  number  of  stations  is  large,  but  as  the  number  of  phases  drops, 
the  3D  model  does  significantly  better.  Further,  because  our  model  is  developed  to  fit  both  P  and  Pn  data,  it 
performs  well  with  regional  data,  teleseismic  data,  or  a  combination  of  the  two. 
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